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Abstract 



We have performed an ecliptic imaging survey of the Kuiper belt with our deepest 
and widest field achieving a limiting flux of m{g')^Q% ~ 26.4, with a sky coverage 
of 3.0 square-degrees. This is the largest coverage of any other Kuiper belt survey 
to this depth. We detect 72 objects, two of which have been previously observed. 
We have improved the Bayesian maximum likelihood fitting technique presented in 
Gladman et al. (1998) to account for calibration and sky density variations and have 
used this to determine the luminosity function of the Kuiper belt. Combining our 
detections with previous surveys, we find the luminosity function is well represented 
by a single power-law with slope a = 0.65 ± 0.05 and an on ecliptic sky density 
of 1 object per square-degree brighter than itir = 23.42 ± 0.13. Assuming constant 
albedos, this slope suggests a differential size-distribution slope of 4.25 ±0.25, which 
is steeper than the Dohnanyi slope of 3.5 expected if the belt is in a state of collisional 
equilibrium. We find no evidence for a roll-over or knee in the luminosity function 
and reject such models brightward of m(R) ~ 24.6. 

Key words: COLLISIONAL PHYSICS KUIPER BELT 



Proposed Running Head: Kuiper Belt Luminosity Function. 

Editorial Correspondence 

Wesley C. Fraser 

Dept. of Physics and Astronomy 

University of Victoria 

Victoria, BC, Canada 

V8W 3P6 



Email address: wesley.fraser@nrc.ca (Wesley C. Fraser). 



2 



1 Introduction 



The study of extrasolar debris and dust disks has revealed that, for at least 
some of these disks to exist as we see them, there must be a source which is 
responsible for replenishment of small-grain dust in the disk. Otherwise, due 
to radiation pressure, the small grain dust would be blown out of their stellar 
systems on time scales shorter than the age of the star. A disk of planetesimals 
embedded in the dust disk which is undergoing collisional evolution is a likely 
source of this dust. Disruptive collisions could produce the necessary influx 
of dust to extend the lifetime of the disk beyond the dust blow-out time. See 
Meyer et al. (2006) for a current review of debris disks. 

The Kuiper belt is analogous to these extra-solar planetesimal disks, and pro- 
vides an excellent laboratory to study and understand the properties of these 
planetesimals and the processes that affect them, including collisional pro- 
cesses, tensile strengths, compositions, and the mechanisms by which they 
formed. Knowledge of the size distribution in the belt can constrain much 
of this information. The size distribution of small objects provides informa- 
tion on the bulk material properties of the objects, (O'Brien and Greenberg, 
2003; Kenyon and Bromley, 2004) while the size distribution of large objects 
can provide information on the conditions under which these bodies formed 
(Kenyon, 2002). 

We performed a 3.0 square-degree survey of the Kuiper belt to determine 
the size-distribution of large (D > 50 km) Kuiper belt objects (KBOs) via 
a measure of the belt's luminosity function. This is the deepest photometric 
Kuiper belt survey of this size ever completed. 

In section 2 we describe our observations. In section 3 we describe our search 
technique and data analysis. In section 4, we derive a relation between the 
size-distribution and the luminosity function. In section 5, we describe the 
statistical analysis used. We present our results in section 6. In section 7 
we discuss the implications of our findings, and in section 8 we present our 
conclusions. 



2 Observations 

The observations used in our survey were taken with the 3.6 m Canada-France- 
Hawaii Telescope (CFHT) and the Cerro-Tololo-Interamerican Observatory 
(CTIO) 4m Blanco telescope. Observations at CFHT were acquired with both 
the CFH12k (0.33 square degree fov.) and MEGAPrime (0.88 square degree 
fov) mosaic cameras while observations at CTIO were made with the Mosaic2 
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camera (0.38 square degree fov), providing 0.67, 0.84, and 1.54 square degrees 
of searchable area respectively for a total of 3.0 square degrees of searchable 
area. Details of the observations are provided in Tables 3 and 4. 

All of these observations were originally acquired for use in searches for satel- 
lites of Uranus and Neptune (Kavelaars et al., 2004; Holman et al., 2004) and 
prior to this work, none of these fields had been searched for KBOs. The ob- 
servations were made when Uranus and Neptune (and all KBOs in each field) 
were at or near opposition, and covered the projected Hill-spheres of each 
planet. The CFH12k and Mosaic2 observations excluded the area closer than 
3 to the planets, to avoid scattered light in the images. These observations 
all occurred near the ecliptic and are well suited to a deep search for Kuiper 
belt objects (KBOs). 

Approximately 10-20 bright, non-saturated stars (~ 20 mag.) per chip were 
used as reference stars for the image reductions. The variation in the aver- 
age reference star magnitudes follow approximately that expected from the 
varying airmass of the observations (see Fig. 1) indicating photometric con- 
ditions during the observations. The remaining scatter is ~ 0.02 — 0.04 mag, 
significantly less than the shot noise for the brightest objects detected. 



3 The Data Processing and Image Search 

To determine the behavior of KBO size distribution requires the discovery of 
a large number of faint moving sources. If long exposures are used, then the 
sources will move far enough (more than the size of the seeing disk) that a 
trailed image results and no additional depth will be achieved. In our deep 
searches we have adopted a strategy of taking exposures short enough that 
trailing losses will be negligible. We then shift the individual exposures to 
account for the expected sky motions of the objects of interest. These shifted 
images are then combined together to achieve the depth needed. To account 
for the various possible sky motions of KBOs, we have shifted the images at a 
variety of rates and angles and then visually examined each of the combined 
images (stacks), searching for point like sources (see Fig. 2). 

3.1 Data Preprocessing 

MEGAPrime: CFHT provides wide field imaging data from the MEGAPrime 
camera in a 'preprocessed' format, ready for science exploitation. The frames 
provided have been processed using the CFHT ELIXIR/FLIPS (Magnier 
and Cuillandre, 2004) processing system. As part of this processing, unique 
bias, flat-field, fringe, and scattered light images are produced for each 
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'camera-run' (typically matching a 14 night dark-run) and applied to all 
frames acquired during that camera run. Dark runs are broken into multi- 
ple camera-runs if a significant change in the camera performance or image 
characteristics is detected. All the MEGAPrime images used in this project 
were acquired within a single camera-run and have been 'detrended' with a 
common set of calibrator images. Standard calibration from CFHT results 
in a flux conserved image with constant sky level across the image that has 
a typical variation of ±3%. 
CFH12K: These data were originally acquired as part of a search for ir- 
regular satellites of Uranus (Kavelaars et al, 2004) and are different from 
the data used for the wide field satellite search presented in Petit et al. 
(2006). Due to the strong time constraints of the imaging, the images were 
acquired in 'classical' observing mode and not automatically detrended by 
the CFHT queued service observing team. In November 2004 the Canadian 
Astronomical Data Center (CADC) acquired the ELIXIR/FLIPS software 
and calibrators for the CFH12K data and subsequently 'detrended' the en- 
tire set of CFH12K images for which global calibrator frames (bias, flats, 
fringe, scattered light) were available. Part of the re-processing effort in- 
cluded the re-processing of the CFH12K frames used in this project. These 
reprocessed images were used in this search and provide a sky flatness of 
typically ±4%. 

CTIO: The CTIO images were originally acquired as part of a project to 
search for satellites of Neptune (Holman et al, 2004). Bias frames were 
acquired on each observing night and a combination of the overscan strip 
and an average of a dozen bias frames was used to remove the instrumental 
ADC bias from each frame. During the period of observations, a number (~ 
15) of independent fields were observed in order to measure the astrometric 
positions of some previously-known Kuiper belt objects. These sequences 
of fields were 'median combined' using the IRAF images. median task with 
high pixel clipping. These combined images provided an excellent flat-field 
frame that was divided into the search images using the IRAF (Tody, 1993) 
mscred.ccdproc task. This process resulted in images which have flat sky 
across the entire mosaic with a typical scatter of ±3% in the sky flux level. 

The curvature of the focal plane for all instruments used in this project was 
small enough compared to the spatial shifts applied to the images during image 
processing, such that spatially flattening the images was unnecessary. 

For the MEGAPrime and CFH12k data, zeropoint calibrations were done 
with the standard Elixir routines and were provided by the CFHT Elixir QSO 
(Magnier and Cuillandre, 2004) and were reported to be Zf ield ^ = 26.46 ±0.02 
mags, and Zfi e id,R = 26.22 ±0.02 mags. Hence these magnitudes are presented 
in the Sloan g' (Smith et al, 2002) and Kron-Cousins R. For the N10032W3, 
N11033, NEP0813NW3, and NEP0815NE3 fields, calibration frames were un- 
available, and a nominal zeropoint for the CTIO mosaic of Z VR = 26.0 was 
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used presenting the magnitudes of the objects discovered in these frames in 
the VR filter. 



3.2 Artificial Object Planting 

To determine the search efficiency as a function of magnitude, artificial objects 
representative of KBOs were added (implanted) in the images. 

The sky rate of motion (9) of an object on a circular orbit observed at oppo- 
sition on the ecliptic, at heliocentric distance r can be approximately given 
by 

arcsec. hr _1 . (1) 

We implanted 100-150 moving sources into each CCD image. All sources were 
added blindly to the data before searching began; the artificial source lists 
were revealed to the operator only after searching was completed. The rates 
and angles of motion of the artificial sources (approximately 1.3 — 4.1 arcsec. 
hr -1 and ±10° from the ecliptic) were typical of objects on circular orbits at 
heliocentric distances from 25 — 100 AU with inclinations between — 70°. 
Each implanted source was given a randomly selected flux, equivalent to that 
of a source between 23 — 27 mag. The distribution of artificial sources is shown 
in Fig. 3. Additionally, five artificial sources with flux levels equivalent to 21st 
magnitude sources were implanted, with zero inclination. These objects have 
sufficient flux to allow us to flag errors in the image combining algorithms 
(failed image subtraction, wrong mask limits, etc.) in advance of searching 
the data, and provided reference moving sources for computing aperture cor- 
rections in the final image combinations. 

To account for image-to-image flux variations due to changes in airmass and 
possible sky transparency, the flux of the planted artificial objects was varied 
with respect to a reference image (usually the middle image in the list) to 
match the average flux variation of the reference stars (see Fig. 1). 

The point spread function (PSF) for each frame was approximated on a per- 
chip basis using the stellar profile of a single bright isolated star. The PSF 
variation across individual images was small enough that the use of a PSF 
that varied with position was unnecessary. 

Artificial sources whose sky motion would result in a drift of the centroid of 
the PSF by more than one pixel during an exposure were represented by a 
series of fainter sources. The number of faint objects was equal to the number 
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of pixels the source would drift during the exposure. The total flux of all the 
fainter sources was equal to the flux of the original source. In this way, we 
fully included trailing effects into the implanted sources, and thus accounted 
for the effects of trailing in our final search efficiencies. 



3.3 Image Subtraction 

During our previous deep searches (Gladman et al, 1998, 2001; Holman et al, 
2004; Kavelaars et al, 2004) we found that the two largest inhibitors to this 
search method are trails in the combined images (caused by bright stars) 
and the enormous human effort required to visually examine the broad range 
of rate/angle combinations needed to ensure detection of KBOs. To combat 
these issues, we utilize an image subtraction routine to remove most stationary 
background sources and implement a new image display method which greatly 
eases the strain on the user during the visual image search. We describe the 
image subtraction routine, and the image display method here. 

A template image was subtracted from each image to remove stationary sources 
from individual images prior to being stacked, thus improving our ability to 
detect moving sources. To create the image subtraction template, an artificial 
skepticism weighted average of all images (Stetson, 1989), with the artificial 
moving sources already added, was created. This method creates a per-pixel 
weighted average, with weights calculated iteratively, and quickly converges 
on an average which places very little weight on spurious pixels such as cosmic 
ray hits. 

For our image subtractions, we chose to use psfmatch3 developed by one of us 
(CJP) for the Supernova Legacy Survey (Pritchet, 2005; Astier et al, 2006). 
This subtraction routine compensates for variations between the image qual- 
ity of the template image and that of the individual images. Several programs 
that perform this task already exist (see Alard (2000) and references therein). 
Psfmatch3 has several advantages. The subtraction kernel can have arbitrary 
form, and does not require representation by a set of basis functions to per- 
form the subtraction. In addition, it can automatically remove both spatial 
and background variations between images. The result of the subtraction rou- 
tine, is smooth subtracted images with zero average backgrounds. A detailed 
description of the method is presented in Appendix A. 

While the image subtraction removed any stationary non-saturated sources, 
cosmic ray spikes and other spurious hot pixels remained. To compensate for 
this, an image mask was applied, such that if a pixel had a value outside a cho- 
sen range, then a 3 pixel by 3 pixel box about that pixel was set to have zero 
value. The lower limit to the range was chosen to be -5 times the background 
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noise (the average background was set to zero in the image subtraction). The 
upper limit was chosen to be ~ 25000 counts, such that most saturated regions 
and cosmic ray spikes were masked out of the data. This procedure masked 
the centres of the brightest KBO sources. This however, did not hinder the de- 
tections of these sources as they were still glaringly obvious even in individual 
frames. 



3.4 Image Stacking 

The subtracted, masked images containing artificial sources were shifted before 
they were stacked, such that each subsequent image was spatially shifted to 
compensate for the predicted motion of KBOs in the frame. Sources whose sky 
motion was well matched by the shifts applied to the images appear nearly 
round in the stack (see Fig. 2) while any residual flux from stationary sources 
was trailed. 

If a source's rate of motion was not well matched by the spatial shifts ap- 
plied to the subtracted images, that source appeared extended in the stack. 
This characteristic trailing provided a very robust means of discriminating be- 
tween real sources and noise. Noise sources were produced during the shifting 
and stacking when positive flux regions from different images were caused to 
overlap by the choice of spatial shifts. These false sources, however, did not 
show the image shape variation characteristic of a real source; false sources 
did exhibit trailing, but of a different length and width compared to a real (or 
artificial) source. As such, false sources were not selected as candidates by the 
trained operator. 

The template image used for the Psfmatch3 subtraction process contains both 
the real and artificial moving sources. These appear as faint trails in the tem- 
plate image. Subtracting this template from the input images results in a low 
flux area behind each moving source. This feature only occurs around the po- 
sitions of real (and artificial) moving sources and provides an additional and 
robust means of source-noise discrimination (see Fig. 2) and was required to 
be present in order to mark a source as a candidate. 

The quality of the image subtraction, and hence the final searchable images 
suffered from a few effects: bad columns in the CCD, bleeding from saturated 
stars, and regions of poor subtraction around bright galaxies. These problems 
were exacerbated by image shifting; the unsearchable area of these regions were 
expanded, reducing the overall searchable area. Gradients of the background 
caused by bad image subtraction around bright sources, produced images that 
were difficult to display, reducing the search efficiency Our efficiency of de- 
tecting artificial, planted sources accounts for all reductions in area as the 
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artificial sources were planted at random locations occupying the full spatial 
extent of each CCD image. 



3.5 Visual Search 



To search the combined stack of shifted images, the stacks were divided up 
into ~ 200 image subsections with a 20 pixel overlap between neighboring 
subsections. One-by-one, a grid of rates and angles like that shown in Fig. 
2 was displayed for each image subsection. Each grid was searched by eye, 
and sources were recognized by their characteristic trailing and subtraction 
wells. We found that a five-by-five grid of shift rates and angles maximized 
the detection efficiency while minimizing the time spent searching. The low 
variations in the sky background, resulting from the template subtraction, and 
an image display tool developed specifically for this searching allowed us to 
rapidly search many data sets. 

The pixel coordinates of potential candidates were recorded along with the 
rate and angle combination that produced the most circular image. A candi- 
date was selected as a planted object if its marked image location was within 
10 pixels of an artificial source location. The list of detected artificial candi- 
dates was used to determine the detection efficiency as a function of artificial 
source brightness. The detection efficiencies for each chip of a given field were 
averaged together to provide a global, per field, detection efficiency which we 
then modeled (see Fig. 4 and Table 5) using the equation 



V (m\r] max ,m^g) 



1 — tanh 



m 



m. 



(2) 



Our deepest and widest field has a rj(m) = 50% threshold at m(g') = 26.4 and 
a sky coverage of 0.84 square degrees. The sky coverage of our combined data 
is 3.0 square-degrees. 

Remaining candidates not marked as artificial source detections were re-examined 
at a finer grid of 8 rates and 8 angles to further discriminate between real 
sources and noise. Candidates were rejected as false positives during this pro- 
cess if the variation of their image shape and trailing length with rate and 
angle did not match that of brighter sources with similar apparent rates of 
motion, whose detection was more robust. While this approach allowed us 
to isolate and remove false detections, the artificial sources were treated dif- 
ferently from the real sources, as we did not examine candidates marked as 
artificial sources in the finer rate and angle grid. The brightnesses of all false 
positives were found to be faintward of the 50% threshold of their discovery 
field, when their fluxes were measured in the same way as all real detections. 
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Therefore, below the 50% threshold, our search efficiency was no longer repre- 
sentative of the true search efficiency, and subsequently, sources fainter than 
the 50% threshold of each field are ignored in our analysis. [ NOTE: Here we 
define the 50% threshold as the point where rj (m) = 50%.] 

A series of follow-up images of the MEGAPrime field as well as the NEP0813NW3 
CTIO field were obtained one or two nights after the initial discovery images. 
Using the short first night discovery arcs, we projected the motion of each 
detection forward to the time of the follow-up observations in order to predict 
the sky location of the source on the second night's images. For 3 of the 25 
sources in the MEGAPrime field, the predicted location placed the sources in 
the gaps between the CCDs of the MEGAPrime MOSAIC. Aside from these 
3 sources, all of our initial detections brighter than the 50% threshold were 
confirmed. Only the faintest initial detection in the MEGAPrime field was not 
confirmed on the second night, likely due to poorer average seeing conditions 
during the second nights observations. All 6 detections in the NEP0813NW3 
field were successfully confirmed on the second night. 

None of our other search fields had follow-up observations. We were confident, 
however, that all detections in the other fields brightward of the 50% threshold 
were real because all detections brightward of the 50% threshold (excluding 
those that fell on chip gaps) in the MEGAPrime and NEP0813NW3 fields 
were confirmed, and all fields were processed and searched with identical pro- 
cedures. 

For the range of rates and angles that we planted objects into the data, we 
found no significant variation in detection efficiency versus rate or angle of 
motion for any of our fields (see Fig. 3). This is important to note because 
any significant variation in efficiency with rate of motion needs to be included 
when deriving the size distribution from the observed luminosity function. 
We did not search for sources with rates of motion consistent with objects 
at distances further than 100 AU as the search method we employed is not 
sensitive to sources beyond this distance; such distant objects would not show 
enough apparent motion over the time-span of our observations to detect the 
faintest sources, and a degradation of the detection efficiency with distance 
would occur. Therefore, our determination of the luminosity function only 
applies to KBOs closer than 100 AU. 



3.6 Source flux measurements and detection confidence 

To measure the flux from each detected source, we stacked the non-subtracted 
images containing the artificial objects. These images were shifted at the rate 
and angle that produced the roundest image for the source whose brightness 
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was to be measured. Three sets of stacks were produced, by averaging the 
images, while using a pixel by pixel cut that threw away pixel values outside 
3-sigma of the mean value. Three averages were made from images from the 
first half, middle half, and last half of the observing period for each field. These 
average images provided three separate, mostly uncorrelated, measurements 
of each source's flux. Aperture photometry was performed using the IRAF 
daophot.phot task, set to use an aperture of 4 pixels radius (close to the size 
of the FWHM of point-sources in our images). Using the five 21st magnitude 
planted 'reference' objects, and the mkapfile IRAF routine, aperture correc- 
tions were determined for each of the three image combinations (first, middle 
last) and were used to correct individual flux measurements for each source. 
The IRAF mkapfile task reported our aperture corrections to an accuracy of 
°"a P er ^ 0.03. The magnitudes, as measured on each exposure set, were aver- 
aged for the final reported magnitudes (see Table 6). We did not measure the 
flux in a particular stack for sources that were within a few seeing disks of 
bright or bloomed stars or galaxies (~ 1 in 10 possible source measurements), 
in that particular stack. 

For two objects in the MEGAPrime data, the measured magnitudes varied 
significantly more than the uncertainty of each individual measurement, indi- 
cating a possible light curve. For these two objects, the magnitudes measured 
from the first image stack were used in our determination of the luminosity 
function, as was done in Bernstein et al. (2004). The second night's data pro- 
vided an additional 3 magnitude measurements of any followed-up source. For 
all but the variable sources, we used the average magnitude measured from 
all image stacks including those from the second night. 

3. 7 Characterization 

As required by the maximum likelihood routine, we parametrize the magnitude 
uncertainty and detection efficiency of our survey. The scatter between each 
artificial source's measured and inserted magnitudes was used to determine 
the precision of our flux measurements. 

For background limited sources, the uncertainty of a source's measured mag- 
nitude can be represented by 

Am = 7 i ( m ^)/ 2 - 5 (3) 

where Z is the telescope zeropoint. 7 is a function of the background b, the 
camera gain g, and the on-source integration time, and is fit to the scatter of 
the artificial source's measured and inserted magnitudes. The fit values of 7 
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for each field can be seen in Table 5. The magnitude measurement scatter for 
each source was then drawn from Eq. 3 using these best fit values. 

The theoretical value of 7 is 



where N is the number of exposures. If we use the typical values from our 
MEGAPrime observations, we expect a value of 7* = 0.12. The measured 
value of 7 = 0.27 is clearly larger. Newberry (1991) however, has shown that 
the noise calculation from Eq. 4 is incomplete, and does not fully account 
for the noise introduced into observations during the data reductions. They 
found that the noise estimate given in Eq. 4 can be too small by a factor of 
2. Similarly, source fluxes were measured off of images that were shifted and 
stacked. If the shifts were slightly different than the source's true motion, the 
measured magnitude would differ slightly from the true magnitude. Thus we 
find that the uncertainties measured from the artificial sources are close to 
the expected values. The net uncertainties used in our luminosity function 
determinations are the 1-sigma shot-noise, Am, aperture correction (a aper < 
0.03), and calibration uncertainties, added in quadrature. 



4 The Luminosity Function 

The size distribution of the KBOs can be determined directly from their lu- 
minosity function. The relation between the luminosity function and the size 
distribution has been derived previously (see, for example, Gladman et al. 
(2001)). We reproduce the relation and discuss complications due to possible 
variations of KBO albedos, using a slightly different approach than that in 
Gladman et al. (2001). 

Consider the number of objects N in a survey field, given by 



where R (r) and S (D) are the radial and size distribution functions, r and 
D are object heliocentric distance and diameter, and A is some normalization 
constant. We model the size distribution as S (D) oc D~ q . The magnitude m 
of an object at heliocentric distance r, geocentric distance A, with diameter 
D and albedo p\ is 




(4) 




(5) 
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m = ^ A + 2.51og 10 (r 2 A 2 D- 2 p^), 



(6) 



where K\ and p\ are wavelength dependent; in R-band, pn ~ 4% gives Kr ~ 
18.8. The smallest observable object at a given distance r in a survey with 
limiting magnitude m max , has diameter D min = p~ 1 ^ 2 r 2 10 ( - Kx ~ Tnrnax ^ 5 where 
we have made the approximation r = A. Similarly, the largest observed object 
at distance r will have diameter D max = p~ 1 / 2 r 2 10^ A_mm4n ^ 5 . 

Using these limits and integrating, Eq. 5 becomes 



, r 2(l-<?) \l0(l-q)(Kx-m min )/5 _ lQ(l-q)(K x -m max )/5 

N(m < m max ) = A R(r)dr ± i - 

J p\ 1 ~ 9 

(7) 

if q 7^ 1. Note: this assumes that D rnin /D max are not the smallest /largest 
objects in the belt; if this were the case, using D min / D max as the integration 
limits as we have defined them here is incorrect. 

For q > 1 and a survey where (q-1) > 2, we have io( 1 ~i)( K - m ^n)/5 ^ 

lQ0--q){K-m max )/5 ^ anc j foe cumulative luminosity function is given by 



N(m < m max ) * A J f^R(r)dr 

* j6d S w^ R ^ dr 1Qa{mmax ~ K) 

where we have substituted q — ha + 1. 

The observed cumulative luminosity function is well represented by a power- 
law of the form 



N{m) = 10 a ("°). (9) 

which gives the number of KBOs per square-degree, on the belt mid-plane, 
brighter than magnitude m. Comparing this with Eq. 8, reveals that the lu- 
minosity function slope a and the size distribution slope q are related by 



q = 5a + l (10) 

From Eqs. 8 and 10 we see that, for a size distribution that is independent 
of distance, the choice in radial distribution is not significant and only affects 
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the interpretation of the normalization of the observed luminosity function, 
not the inferred size distribution slope. 

Here we have assumed that the KBO albedos do not vary with distance, or 
size. If the distribution of KBO albedos varies only with distance, then this 
will only affect the normalization of the luminosity function. 

The albedo of Pluto and Eris (D > 2000km) are ppi uto ~ 0.6 and Pehs ~ 
0.6 - 0.86 (Young et al, 2001; Bertoldi et al, 2006; Brown et al, 2006) while 
smaller objects (D ~ 100 km) have been shown to have albedos p ~ 0.06 
(Grundy et al, 2005). These data suggest an increase in albedo with size. We 
can understand the effects of such a trend on the interpretation of the LF by 
considering a toy model where KBO albedos vary as p ~ D~P. This functional 
form retains the analyticity of Eq. 5, and reveals the effects of albedo variations 
on the inferred size distribution. Under this assumption, Eq. 8 becomes 



4(1-9) 



N(m < m max ) = A J R(r)drl0 . (11) 



Q 

Thus the slope of the size distribution is 



Q = 5« (1 - f J + 1 (12) 

We see that in the case of constant albedo (j3 = 0) we get back Eq. 10. 

Current albedo data imply /3 ~ — 1 (Stansberry et al, 2007). Thus, an estimate 
of q which assumes (5 — potentially under-estimates the steepness of the 
intrinsic size distribution. Our knowledge of the relation between albedo and 
size however, is currently insufficient to constrain (5. 

The reader is cautioned that the current determinations of the size distribution 
from the Kuiper belt luminosity function are based on a few poorly constrained 
estimates of KBO physical properties and the inferred slope q is probably 
underestimated. 



5 Maximum Likelihood fits 



To determine the optimal values of a and m , we use a maximum likelihood 
method. We extend the maximum likelihood analysis of Gladman et al. (2001) 
to account for observations made in different wave-bands, as well as systematic 
differences in the normalization, m , that may affect the results of combining 
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separate data sets, such as systematic calibration errors, sky density variations, 
and variations in average colour of the observed KBOs compared to the average 
colour of all KBOs. 

Given that the probability of detecting an individual KBO is independent of 
detecting the next, the likelihood function for a single survey k takes the form 



L k (a,m \m 1 ,m 2 , ...) oc exp Nk J|Pj (13) 

i 

where mj is the magnitude of detection i, N k is the number of objects expected 
to be detected in a given survey, and p is the probability of having object i 
given the underlying luminosity function. N k is given by 

N k = J dm Vt rj(m) E(m|a, ra ) (14) 

where Q is the survey area, r)(m) = rj(m\rj max , to*, g) is the detection efficiency 
for an object with magnitude to, and S(to|q;,to ) is the differential density of 
objects on the sky. p and E(to|q;,to ) are given by 

= £(,»!«, ,»„),,(,»). (15) 

and 



E(m\a, m ) = ^Sllll = l n (i ) a 10^"°), (16) 
dm 

where a is the slope of the luminosity function (LF), m Q is the magnitude at 
which the sky density is 1 KBO per square degree, and €j(m) is a functional 
representation of the photometric uncertainty for object %. A formal derivation 
of this likelihood function is presented in Loredo (2004). 

We choose to represent the magnitude uncertainties as gaussian. ie. 



q(to) = e^? (17) 

^ttAtoI 

where Atoj is the uncertainty in the magnitude measurement of each object. 
This treatment of uncertainties for faint sources is incorrect, but is a suffi- 
cient approximation for detections brightward of the 50% efficiency threshold. 
Brightward of this threshold, the gaussian approximation will not affect the 
results of the maximum likelihood inference (Bernstein et al, 2004). 
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For a group of surveys with calibration uncertainties and colour offset vari- 
ations much smaller than the uncertainty in the flux measurements of the 
brightest objects detected, the net likelihood resulting from combining multi- 
ple surveys together is the product of the likelihoods of each individual survey. 
In reality, photometric calibration and colour offset uncertainties are not in- 
significant, and therefore, need to be considered when combining different 
surveys. 

Additionally, the apparent sky density of KBOs and, hence, m are strong 
functions of latitude and longitude (Kavelaars et al, In Press). When com- 
bining separate surveys, differences in the flux-limited sky densities will skew 
the inferred slope of the LF. 

We account for these effects with the addition of two new parameters for each 
survey, the colour parameter C k , and the density parameter Am Qtk . Eqs. 14, 
15, and 16 become 



N k = J dm Q r](m\r} max ,m*,g) S fe (m - C k \a,m , Am 0)fc ), (18) 
P i = J dmT, k (m- C k \a,m , Am 0jk ) e^ra), (19) 



and 



£(m - C k \a, m , Am Qyk ) = ln(10) a io«(™- c *-(™°- A ™°.*)). (20) 

As can be seen, for the case of a power-law LF, the parameters C k and Am 0tk 
are degenerate. For a non power-law model however, C k and Am 0tk are no 
longer degenerate, and need to be treated separately. As we consider non 
power-law LFs in this manuscript, we choose to treat C k and Am 0tk as different 
parameters when evaluating the likelihood. 

This treatment substantially increases the number of parameters in the fit, 
while not providing any new information about the shape of the LF. Thus 
we treat these additional parameters as nuisance parameters, and marginalize 
(integrate) the likelihood equation over the expected ranges of each parameter. 
The final likelihood equation when combining M surveys is 



at 

L(a,m ,C 1 ,C 2 , ...,CV, Am 0j i, Am 0i2 , Am ^) = JJ / / L k (a,m ,C k , Am 0ik )dC k dAm 0tk . 

k=l 

(21) 



16 



6 Results 



Each detection in our survey is listed in Table 6 along with estimated barycen- 
tric distance and inclination determined from fit_radec (Bernstein and Khusha- 
lani, 2000). The fit_radec routine determines a set of orbital parameters that 
best fit the observations in a least-squares sense. Because of the large degen- 
eracies between orbital parameters for short-arc observations such as those 
presented here, fit_radec initially assumes that the objects are near perihelion, 
and are on bound, nearly circular orbits, near the ecliptic plane. The routine 
determines 2-sigma uncertainties, which contain 95% of the orbits which are 
consistent with the measured positions of the objects. In this survey, we have 
detected 72 KBOs, 53 of which are brightward of the 50% threshold of our 
search, and these 53 detections are used in our maximum likelihood analysis. 

As the UNE and UNW fields are close on the sky, observed together through 
the same filter/telescope, and the skies were photometric during those obser- 
vations, we combine the two fields for the maximum likelihood fits. Similarly, 
we treat the N10033 and N10032W3 and the NEP0815NW3 and NEP0815NE3 
fields in the same fashion. We refer to these combined fields as the UN, 
CTIO01, and CTIO02 fields. The CTIO01 and CTIO02 fields were not com- 
bined together, as they were observed in separate years. 

To extend our results, we include the observations from other surveys that 
have well-measured efficiency functions for each field in the survey. We define 
the F07 sample as those objects detected in the surveys from Jewitt et al. 
(1998), Gladman et al. (1998), Chiang and Brown (1999), Gladman et al. 
(2001), Allen et al. (2002), and Petit et al. (2006), as well as those detected 
in the UN, CTIO01, and CTIO02 fields. 

In the F07 sample we also include all on ecliptic detections in Trujillo et al. 
(2001) (T01) in fields with detection efficiencies classified as 'good' or 'medium'. 
Because the T01 images cover a large range of latitudes, longitudes, and detec- 
tion efficiencies, we divided the survey into smaller "sub-fields" to avoid large 
density variations from field to field. Each sub-field spans ~1 - 1.5 hrs. RA. We 
do not include any of the high-latitude data from T01 as the KBO latitude 
distribution is not sufficiently understood to be able to predict the density 
of KBOs at high latitudes and we are thus unable to provide a reasonable 
estimate of Am Dife for the high-latitude fields. We avoid complications due to 
detection efficiency variations by further sub-dividing the data from T01 into 
separate groups based on the good and medium detection efficiencies defined 
in that work. We do not consider any T01 fields with detection efficiencies 
classified as 'poor'. The sub- field divisions listed by field ID (see Table 2 in 
Trujillo et al. 2001) are presented in Appendix B. 
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We include all fields from the above surveys, even those with few or no detec- 
tions in the F07 sample. This is the sample on which all our LF analysis was 
performed. 

We chose to exclude the data from Bernstein et al. (2004) as the intent of this 
work is to determine the shape of the bright end of the luminosity function. 

When computing our LF estimate, we ignored all detections fainter than the 
77(771) = 50% detection threshold of each individual survey, and truncated 
( "cut" ) the detection efficiency to zero faintward of that threshold. Correctly 
determining the efficiency function at low levels is notoriously difficult and 
prone to error. Effects such as Malmquist-bias further complicate detection 
efficiency measurements, at low brightnesses. 

To test for a bias in the fits caused by an efficiency truncation, we ran a series 
of Monte-Carlo realizations that simulated the observation of a luminosity 
function from a single survey. We also performed this analysis using a set 
of surveys with parameters like those in the F07 sample. We performed 2000 
random realizations for efficiency cuts ranging from 0-90%. We found that, for 
a single survey, the mean of the estimated LF slope is biased to steeper slopes 
when truncating the detection efficiencies at some non-zero threshold. The 
bias was found to decrease asymptotically to zero as the efficiency threshold 
was moved faintward. 

We found the bias caused by including an efficiency cut is removed however, 
when combining multiple surveys if the separate surveys have different detec- 
tion efficiencies. Most detections in a survey occur where 77 (to) ~ 80%. Thus, 
when two or more surveys are well separated in their limiting magnitudes, 
so are their detections, which creates a "lever arm" that dominates the slope 
determination, effectively removing any efficiency cut bias in the fit. 

For each field with 5 or more detections which is the minimum number of de- 
tections required to provide a reliable measure of the LF, we determined the 
best-fit parameters for the LF using Eq. 13 (see Table 7). The 1, 2, and 3-sigma 
credible regions of these fits are presented in Fig. 5. The 1-sigma uncertain- 
ties of the best-fit parameters, taken as the extrema of the 1-sigma credible 
regions, are presented in Table 7. Presenting the fit parameter uncertainties in 
this way describes the full range of allowed values for each parameter individ- 
ually. These uncertainties are somewhat misleading as they do not describe 
the correlation between a and m (see Fig. 5). 

The average (V-R) colour of KBOs is < V — R >= 0.56 (Hainaut and Delsanti, 
2002). Using the VR and R magnitude transformation given by Allen et al. 

(2001) , as well as the relations between V, R, r' and g' given by Smith et al. 

(2002) , we determined the expected offsets between the various filters used in 
the surveys considered here. We found that the average colours of KBOs are 
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<VR-R>= 0.03, < r' - R >= 0.26, and < g' — R >= 0.95. Using these 
colours to offset the separate fields to R-band, we find the estimates of m 
provided independently by each survey are not consistent with a single value 
(see Fig. 5). Small shifts in m are necessary to make all fields agree at the 
1-2 sigma level. This is acceptable when considering possible sky density 
variations, and photometric calibration errors (see Section 5) and justifies the 
more complicated likelihood equation given by Eq. 21. We therefore shifted 
all data to R-band using the expected average KBO colours, and used Eq. 21 
when determining the best-fit LF parameters. 

To determine the possible range in m values between surveys, we consider 
a toy model of the Kuiper belt. As the majority of the objects observed in 
the surveys considered in this work are Plutinos and classical belt objects 
(CKBOs), we only consider these objects in the model. We represent the LFs 
of both populations as power-laws with the same slope, but with different 
observed sky densities. Then the observed cumulative LF (assuming constant 
detection efficiencies) is a power-law, and is given by 

jV {, = lQ a ( m ~ m °^) _|_ JXQ a ( m ~ m o,p) = ^ga(m-m„) ^2) 

where m 0jC , m 0iP , and m are the magnitudes at which one object per square 
degree is observed for the CKBOs, the Plutinos, and all populations respec- 
tively and / is the ratio of the number Plutinos to CKBOs. 

The mean eccentricity of Plutinos is e = 0.15, and they are ~ 30% as numerous 
as the classical belt objects of the intrinsic population, ie. / = 0.3 (Kavelaars 
et al. (In press)). The magnitude of an object as a function of heliocentric 
distance r is given by m ~ K + 10.01og(r). Thus the variation in m 0tP for 
Plutinos at perihelion versus aphelion is ~ 1 mag. If all CKBOs are observed 
at r = 43 AU, then the difference in m between fields that observe Plutinos 
at perihelion versus aphelion is ~ 0.8 mags. This is an upper estimate of the 
true offset because not all Plutinos come to perihelion at once, but provides 
a reasonable range of integration for each m j, parameter. 

Latitude differences between surveys will cause variations in m requiring 
knowledge of the latitude distribution of the Kuiper belt. While the location 
of the mid-plane is unclear, it is secure that the mid-plane is inconsistent with 
the ecliptic (Brown and Pan, 2004; Elliot et al, 2005); there is disagreement 
about whether the mid-plane is the invariable plane. Brown (2001) has shown 
that the KBO latitude distribution above the plane has a Full- Width-Half- 
Maximum of ~ 7°. Most of the observations considered in the F07 sample are 
made near the ecliptic. If the true plane of the belt is similar to the invariable 
plane, or the Kuiper plane proposed by Brown and Pan (2004), then the den- 
sity variation due to latitude between ecliptic fields is at most 5% and would 
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affect the m ^ values by ~ 0.05 mags. This variation is small compared to 
that expected from changes in longitude, and is ignored in our analysis. We 
bound the integration of Eq. 21 over the m ^ parameters between ±0.4 mags. 

The standard deviation of the < V — R > colour from the MBOSS data set is 
o~(V — R) ~ 0.15 (Hainaut and Delsanti, 2002). Hence, the expected variation 
in average object colour between separate surveys due solely to KBO colour 
variation is ~ 0.15/ ^/N where N is the number of objects in a survey. For 
the fields considered in the F07 sample, N ~ 1 — 15. Zeropoint calibration 
errors between separate surveys have been as high as 0.1 magnitudes. Thus 
a reasonable range of offsets due to calibrations and colour variation is ~ 
0.2 mags. This offset can occur in either a positive or negative sense on the 
measured magnitude. Thus we bound the integration of Eq. 21 over the Ck 
parameters between ±0.2 mags. 

Maximizing Eq. 21 using the F07 sample, we find a best fit slope of a — 
0.65 ± 0.05 and normalization m = 23.42 ± 0.13. The likelihood contours of 
the fit are shown in Fig. 6. 

The factor e~ N in Eq. 21 weights the fit towards the lowest possible number of 
detections. The m that maximizes the likelihood equation given by Eq. 21 is 
the maximal value within the range allowed by the Am 0i k that best describes 
each individual field considered. Thus, the best-fit m Q is not applicable to any 
one field, but rather is a value typical of all data considered in the fit. 

To test if the fit is consistent with the observations, we employed a series 
of Monte-Carlo simulations. The simulations involved random realizations of 
a number of objects drawn from the best-fit power-law model (a = 0.65, 
m = 23.42) equal to the number of objects in the F07 sample with the ran- 
dom magnitudes scattered according to our uncertainty model (see Eqs. 3 and 
17). A best-fit power- law of each realization was determined using our maxi- 
mum likelihood technique, and the distribution of maximum likelihoods was 
determined from 1000 realizations. We found that the probability of getting a 
maximum likelihood less than or equal to the maximum likelihood computed 
from the F07 sample was 43.2%: the model is fully consistent with these data. 

In Fig. 7 (top) we present a histogram of the differential LF of the F07 sample 
corrected for detection efficiencies. These data are well described by our best 
fit. The reader is cautioned however, from drawing conclusions about the LF 
from this diagram alone. All source magnitudes have been adjusted to R-band 
using typical KBO colours (see above). The observational data from different 
surveys however, contain calibration and colour offsets, and the observations 
have not been adjusted to reflect variations in sky densities, as a standard 
model of the density variations is not known. The minor discrepancies between 
the observed LF and the model apparent in Fig. 7 maybe caused by these 
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effects. 



In Fig. 7 (middle) we show the net effective area for all fields as well as a 
differential histogram of the object magnitudes shifted to R-band with com- 
pleteness corrections. The different fields which are included in the F07 sample 
all have different 50% thresholds, resulting in the broad fall off in effective area 
between m(R) ~ 21 to 26. 

In Fig. 7 (bottom) we show the logarithm of the ratio of the observed and fit 
differential LFs, as a function of magnitude. This figure shows that the LF 
from m(R) = 21 to 26 is well described by a single power-law. 



6.1 Broken Power-law Models 

Bernstein et al. (2004) conclude that the slope of the luminosity function 
rolls over to shallower values for fainter magnitudes and that the luminosity 
function is well described by a rolling power-law given by the functional form 

£(m) = S23 io Q ( m - 23 )+«'( m - 23 ) 2 (23) 

where S(m) is the differential surface density of objects, S 23 is the surface 
density of objects at m(R) = 23, a is the bright end slope, and a' is the slope 
derivative. They found that a rolling power-law was a better fit to their ob- 
servations than the single power-law used here. This suggests that a deviation 
in the form of a flattening of the power-law might be visible in our data. 

To look for evidence of a roll-over in the KBO LF we fit Eq. 23 to the F07 
sample. As the density parameter for Eq. 23, S 2 3, is different from m in 
Eq. 9, we introduce a density offset parameter as a multiplicative factor 
in front of Eq. 23. We maintain the same range in colour and density offsets 
as used in the power-law fit (see above) and marginalize and AS fc over the 
range ±0.2 and ±0.25 respectively. Our maximum likelihood method gives a 
best-fit A = 0.79 ± 0.14, a = 0.74 ± 0.09, and a' = -0.03 ± 0.04. Bernstein 
determined a best-fit of A = 1.07, a = 0.66±0.03 a' = -0.05±0.015 (1-sigma 
uncertainties were extracted from Fig. 4 of Bernstein et al. (2004)), consistent 
with our results. The increase in uncertainty of a and a' in our fit versus that 
from Bernstein et al. (2004) is caused by the marginalization in Eq. 21. We feel 
that our approach provides a more realistic estimate of the true uncertainties. 

Additionally we fit a simple broken power-law model by Eq. 21. The model 
LF has a sudden change from the bright-side slope ct\ to the faint-side slope 
a 2 at a break magnitude m B , and is given by 
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if m < tub, 
if m > iris. 



(24) 



The best-fit parameters from maximizing the likelihood for the broken power- 
law LF while marginalizing over nuisance parameters Ck and Am 0t k as per 
Eq. 21 are m = 23.2 ± 0.5, a x = 0.69 ± 0.08, a 2 = 0.57 ± 0.2, and m B = 
24.4 ±0.7. 

The improvement of the maximum likelihood value when considering the 
rolling and broken power-law LFs given by Eqs. 23 and 24 fits is only a few 
percent over that of the single slope model. Hence the additional degrees of 
freedom included in the more complicated LFs compared to a single slope 
power-law does not substantially improve the fit and the higher order mod- 
els are not statistically warranted in the range of magnitudes considered in 
the F07 sample. Indeed the best fits of both equations are consistent at the 
1-sigma level with no break at all (a' = and a 2 = «i). 

To test at what magnitude a roll over to shallower slopes would be incon- 
sistent with the F07 sample, we made use of Monte-Carlo simulations and 
the Kolmogorov-Smirnov (KS) test. In these simulations, we simulated the 
observation of a broken power-law LF of the form of Eq. 24, with cci = 0.65 
breaking to some faint-end slope a 2 at break magnitude m B . The simulated 
surveys and total number of detections were chosen to match the surveys and 
number of detections in the F07 sample. We included the effects of variable 
m and colour offsets between surveys by randomly selecting these offsets for 
each individual survey. 

For a given break magnitude and faint-end slope, a parent population of 
~ 10, 000 objects was generated. This parent population was simulated with 
a set of random Ck and Am o fc offsets linearly sampled from the marginaliza- 
tion range described above (±0.2 for and ±0.4 for Am 0i (.). Calibration and 
colour offsets occur in the observations in the F07 sample. The random sam- 
pling we implement provides a simple means of generating offsets consistent 
with the way in which those observations were made. 

From the parent population, a sub-sample of objects was bootstrapped and 
the KS-statistic of this sub-sample compared to the parent population was 
calculated. This procedure was repeated for 1000 random sub-samples of the 
parent population. In this way, the KS statistic distribution was bootstrapped 
from the parent population. 

The KS statistic of the F07 sample when compared to the parent population 
was compared to the distribution of KS-statistics generated from the boot- 
strapped samples. The roll-over model (a 2) m B ) was rejected if 99.7% of the 
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bootstrapped KS statistics were smaller than the KS statistic of the F07 sam- 
ple. This was repeated 25 times using different random realizations of the Ck 
and Am 0) jt offsets for each choice of faint-end slope and break magnitude. The 
repetition was necessitated by the randomness included with variable Ck and 
Am 0j fc values; in certain circumstances, a particular set of colour and m offsets 
will enhance the chances of observing a break. A model could not be rejected 
if the break was not rejected in 2 or more of the 25 simulation repetitions for 
that model. 

Presented in Fig. 8 is the 95% probability contour that a particular break 
model is rejected by the F07 sample. The contour has been smoothed to remove 
the effects of course sampling in a 2 -' m B space. As expected, brightward of a 
particular magnitude most break models would likely have been detected in 
the observations. From these simulations, we conclude that the observations 
cannot reject the possibility of a break in the Kuiper belt LF with break 
magnitudes fainter than m B (R) ~ 24.3. 

Bernstein et al. (2004) find that the LF is well described by the harmonic 
mean of a steep power-law for bright objects and a shallow power-law for faint 
objects with both power-laws contributing equally at magnitude R eq . They 
find R eq ~ 22.8 — 23.6. This is not however inconsistent with the results of 
our KS test as and R eq are not equivalent parameters between the two LF 
models. The reader is cautioned about drawing conclusions from comparison 
of these two parameters. 

While no break is apparent in the observations, the best-fit single sloped model 
is inconsistent with the results from Bernstein et al. (2004) at m(R) > 28 as 
the number of objects detected in that survey are too few by a factor of ~ 6 
from that predicted by the model. To account for this, the LF must roll-over 
to shallower slopes at some faint magnitude not present in the F07 sample. 



6.2 The Size Distribution 



Our best-fit luminosity function slope is a = 0.65±0.05. Under the assumption 
that Eq. 10 holds, we find that the KBO differential size distribution has a 
slope of q = 4.25 with a 1-sigma uncertainty of 0.25. If the belt is in a state 
of collisional equilibrium, we would expect the slope to be q ~ 3.5. This 
is inconsistent with the inferred size distribution at the 3-sigma level. We 
conclude that, for objects larger than r pa 50 km, the size-distribution of the 
belt is inconsistent with a system in collisional equilibrium. 
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7 Discussion 



To understand what the KBO size distribution tells us about the history of 
objects in the outer solar system, we must interpret our observations in terms 
of models that account for the size distribution of a belt of planetesimals. Gen- 
erally, there are two broad types of models that attempt this, fragmentation 
models and accretion models. 

Analytic fragmentation models are those in which a series of equations ac- 
counting for accretion and collisional disruption are solved, producing a steady- 
state collisional-cascade equilibrium. These models vary in complexity. Some 
assume that each disruption produces a number of equal size collision remnants 
(Dohnanyi, 1969; Pan and Sari, 2005), while some model object-disruptions 
with a distribution of collision remnant sizes (O'Brien and Greenberg, 2003). 
Some calculations also include various models of KBO physical strengths 
(O'Brien and Greenberg, 2003). The detailed results depend on the calcu- 
lations, but generally the outcome is a differential size distribution in a quasi- 
equilibrium steady-state with a slope equal to or less than Dohnanyi slope, 
q < 3-5. 

Pan and Sari (2005), assumed that the strengths of large Kuiper belt objects 
are gravity dominated and modeled a collisional cascade in which all collisions 
were purely disruptive (no accretion). They found that, for objects smaller 
than some break size, the size distribution slope was q = 3 and could not 
account for the steep slope observed for large KBOs with this model. 

O'Brien and Greenberg (2003) accounted for accretion and fragmentation, and 
a variation of object strength as a function of size. They found that the equi- 
librium size distribution slope was related to the power-law slope describing 
the variation in object disruption energy per unit mass with size. From this 
model, the Kuiper belt size distribution inferred from the F07 sample implies 
a unit disruption energy that decreases rapidly with increasing size. This is 
consistent with objects whose physical strength is dominated by tensile forces, 
but is inconsistent with objects whose strength is dominated by gravity; these 
objects have an increasing disruption energy with size. Objects in the tensile 
strength dominated regime are typically smaller than 1 km in size. The objects 
we observe here are likely to be gravity dominated and the scaling of KBO 
disruption energy versus size implied by the model from O'Brien and Green- 
berg (2003) seems unlikely. We conclude from these results that the observed 
LF and the inferred size distribution is inconsistent with analytic equilibrium 
models. 

The large diversity of fragmentation models have similar outcomes; they pre- 
dict a shallow size distribution slope of q ~ 3.5 which is incompatible with 
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the large object size distribution inferred from the luminosity function of the 
Kuiper belt. 

Numerical accretion models are those in which the size distribution, and orbits 
of a population of bodies is calculated. Unlike analytic collisional cascade 
models, the size distribution calculated from accretion models is not assumed 
to be in steady state, but evolves in time along with the orbital distribution. 
Models which account for accretion and fragmentation of planetesimals in the 
region of the Kuiper belt, predict a broken power-law size distribution with a 
steep slope for large objects (D > 10 km) that rolls over to a shallower slope 
for smaller objects (D ~ 2 km) (Kenyon and Bromley, 2001; Kenyon, 2002; 
Kenyon and Bromley, 2004). 

These models have two general evolutionary phases. The first phase is plan- 
etesimal accretion, in which a large reservoir of planetesimals on nearly circular 
orbits accrete to form larger bodies via low encounter-velocity collisions. This 
process rapidly produces a steep size distribution for large objects (q > 5). 
Models from Kenyon and Bromley (2001) and Kenyon (2002) calculate planet 
growth in the 40-47 AU zone. These models have a relatively low-mass intial 
Kuiper belt, and do not include stirring from Neptune. They calculate planet 
growth before Neptune attains its current mass and orbit, and start from bod- 
ies with radius < 1km. Initially, the size distribution of large objects is very 
steep (q > 5). After ~ 100 — 300 Myr, the slope of the largest objects flattens 
to q ~ 4.1 - 4.4. 

The second phase is started when the belt members are stirred up onto dynam- 
ically excited orbits (via interactions with Neptune, oligarchs, etc.) such that 
mutual collisions are mainly catastrophic. Collisional grinding rapidly reduced 
the slope of the size distribution for small bodies (D < 1 km). The longer the 
collisional evolution continues, the larger the radius at which the break to 
shallower slopes occurs. In the models from Kenyon and Bromley (2001) and 
Kenyon (2002) the size distribution breaks to q < 3.5 at D ~ 2 — 20km after 
~ 1 Gyr. 

More complicated models that start with a more massive initial Kuiper belt, 
and include effects of gas-drag and a more detailed model of object strength, 
produce large object slopes of q ~ 3 and after 4.5 Gyr, and a break at D ~ 
2 — 40 km to very flat slopes of q ~ — 0.5 for small objects. (Kenyon and 
Bromley, 2004) 

The observed size distribution exhibits a steeper slope than those determined 
by Kenyon and Bromley (2004) and is more consistent with models that start 
with a less massive Kuiper belt. In these models, the break radius grows 
larger and the large object slope flattens as accretion and collisional evolution 
continues. An accurate measure of the break size in addition to the large 
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object slope determined here, combined with modeling, may further constrain 
the duration of accretion in the region of the Kuiper belt. 

Accretion models which include migration of Neptune can account for some 
of the features of the KBO orbit distribution, but come at the cost of com- 
plicating the accretion process (Kenyon et al. 2007 and references there in). 
In these models, most of the current Kuiper belt consists of objects that are 
scattered by Neptune to their current positions, and originate from regions 
interior to Neptune's current orbit. Accretion of the KBOs in their original 
locations proceeds until stirring and scattering by Neptune disrupts the ac- 
cretion process. 

In these migration models the KBOs are initially closer to the Sun than there 
current locations and accrete in regions of high density, and therefore both the 
accretion and collision process occur on much faster time-scales than in-situ 
formation models predict. The migration time-scale becomes an important 
parameter in these models; too fast a migration and the largest KBOs are not 
formed, while too slow, and the break size can evolve to be large enough such 
that it would be detected in the available observations. 

Numerical accretion models of KBO formation, are generally in better agree- 
ment with the observed large size distribution than collisional cascade equilib- 
rium models. Accretion models suggest that knowledge of the size distribution 
slope, and the radius at which the break to shallower slopes occurs will provide 
strong constraint on the accretion dynamics, and KBO formation time-scales. 
Due to the large degree of complexity in these models however, there is much 
work needed to interpret the observed size distribution. Models which achieve 
collisional-cascade equilibrium for large objects cannot account for the ob- 
served steep size distribution. 



8 Conclusions 

We have performed a survey of the Kuiper belt with an sky coverage of 3.0 
square degrees with our deepest field having a depth of m(g') = 26.43 magni- 
tudes. An analysis of current survey data confirms that the luminosity function 
of the belt is well described by a single slope power-law between m(R) = 21—26 
with a slope of a = 0.65 ± 0.05 and normalization m = 23.42 ± 0.13 which is 
typical of all the fields considered in the F07 sample. 

We have shown the necessity of considering calibration, and density effects 
when inferring the luminosity function from the combination of different sur- 
veys that have observed separate regions of the sky, and are not directly cal- 
ibrated to one-another. Such effects can skew the inferred LF away from the 
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true form, and cause uncertainties in the inference to be highly underesti- 
mated. 



We conclude that: 

(1) The Kuiper belt is not in a state of collisional equilibrium for objects 
larger than D 50 km. 

(2) There is no evidence for a change in the slope of the luminosity function 
for magnitudes brighter than m(R) ~ 24.3 corresponding to objects with 
diameters D > 110 km at 40 AU. The sample of observations considered 
here is consistent with a best-fit model that breaks at m(R) ~ 24.4 to a 
slope of a 2 ~ 0.6, but the more complicated fit is not warranted by the 
slight improvement in the maximum likelihood value. 

(3) The observed slope of the luminosity function implies that the size dis- 
tribution of the Kuiper belt is consistent with a power-law with slope 
q = 4.25 ± 0.25 for objects with diameters larger than D ~ 50km. 
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A Psfmatch3 

Psfmatch3 is a subtraction routine that compensates for variations be- 
tween the image quality of the template image and that of individual 
images. A description of the basic algorithm follows. 

Consider an observed image (i, j refer to pixel row and column), and 
a suitably chosen reference image 1^. We define an x convolution 
kernel Kij such that the convolution K * I provides some "best" (in a 
least-squares sense) estimator, O, of O; that is, the optimal kernel K is 
the one that minimizes Yri,j(Oij — Oij) 2 . Differentiating with respect to 
each kernel element yields a system of M = n\ linear equations 
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which can be solved for the n\ kernel elements. 

This basic procedure can be improved considerably by adding an extra 
term for sky differences, viz. 



O = K*I + As, 

where As is a constant denoting the difference in sky between the two 

V\ . (Oij-dij) 2 

images. Weighting by errors is easily included by minimizing — — — T . 

Perhaps most important is allowing for spatial variations in the kernel 
(and sky background); this is accomplished simply by solving for poly- 
nomial coefficients representing the spatial variation of each kernel coef- 
ficient. 

There are a number of features and advantages to this method of 
matching the point spread function on two images. 

• The kernel K has arbitrary form; it does not need to be symmetric, and 
can handle PSF variations that may not always be handled by methods 
involving a basis set of functions to represent the kernel. The kernel 
can even be solved to perform deconvolution (which is the case when 
the reference image is erroneously chosen to have worse seeing), though 
this results in noise amplification. 

• K is not necessarily normalized to unity; this automatically takes out 
transparency fluctuations. 

• The method automatically removes small shifts (even spatially variable 
shifts) between two images. 

• As noted above, the method can easily be adapted to include back- 
ground differences, spatially variable backgrounds, and spatially vari- 
able kernels. 
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2 Field Divisions from Trujillo et al. (2001) 



Table 1 
Field Details. 



" 11 

Field 


R.A. (hrs.) 


T71 1' A.' T A.' A. 1 1 (~)\ 

Ecliptic Latitude (°) 


Area (Sq. ) 


# Objects 


t~i m 

Emciency 


TE1G 


8-9 


-0.5 - 0.5 


7.59 


15 


Good 


mn — \ -1 TV If 

TE1M 


8-9 


-0.5 - 0.5 


3.3 


4 


Medium 


TE2G 


10 - 11 





4.35 


11 


Good 


TE3G 


11 - 12.5 





6.27 


12 


Good 


TE3M 


11 - 12.5 





2.31 


4 


Medium 


TE4G 


~21-23 





2.64 


4 


Good 


TE4M 


~21-23 





2.97 


3 


Medium 


TE5G 


~23 - 1 





0.66 





Good 


TE5M 


~23 - 1 





4.39 


7 


Medium 
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Table 2 

Trujillo and Brown (2001) pointings included in each field division. 



Field 


Pointings 


TE1G 


476727o, 476728o, 476729o, 476848o, 476849o, 476850o, 476851o, 476852o, 476853o, 
476854o, 476855o, 476856o, 476857o, 476858o, 476859o, 476984o, 476985o, 
476986o, 476987o, 476988o, 476989o, 476990o, 476991o 


TE1M 


476717o, 476718o, 476719o, 476720o, 476721o, 476992o, 476993o, 476994o, 476995o, 
476996o 


TE2G 


476885o, 476886o, 476887o, 476888o, 476889o, 476890o, 476891o, 476892o, 476893o, 
476894o, 476895o, 476896o,476924o 


TE3G 


476758o, 476759o, 476760o, 476761o, 476762o, 476763o, 476764o, 476765o, 476766o, 
476767o, 476768o, 476769o,476795o, 476796o, 476797o, 527174o, 527305o, 
527458o, 527461o 


TE3M 


476798o, 476799o, 527175o, 527306o, 527307o, 527459o, 527460o 


TE4G 


502047o, 502048o, 502049o, 502183o, 502184o, 502215o, 502217o, 502218o 


TE4M 


502050o, 502051o, 502052o, 502182o, 502185o, 502186o, 502214o, 502216o, 502374o 


TE5G 


502102o, 502139o 


TE5M 


502098O, 502099o, 502100o, 502101o, 502103o, 502136o, 502137o, 502138o, 502140o, 
502248o, 502249o, 502250o, 502426o 
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Table 3 
Field Details 



Field 


Telescope 


Date (UT) 


Camera 


Filter 


Exp. Time 


a 


5 


Seeing (") 


UNE 


CFHT 


2001-08-24 


CFH12k 


Kron-Cousins R 


16 x 480 s 


21:41:06.6 


-14:28:04.9 


1.0 


UNW 


CFHT 


2001-08-25 


CFH12k 


Kron-Cousins R 


17 x 480 s 


21:38:34.9 


-14:37:15.6 


0.87 


MEGA 


CFHT 


2004-09-15 


MEGAPrime 


g' 


48 x 240 s 


22:24:46.4 


-10:46:51.3 


0.74 




CFHT 


2004-09-16 


MEGAPrime 


g' 


45 x 240 s 


22:24:37.9 


-10:47:39.0 


0.83 


N10032W3 


Blanco 


2001-08-10 


Mosaic2 


VR 


40 x 480 s 


20:38:35.8 


-18:01:02.5 


1.4 


N10033 


Blanco 


2001-08-11 


Mosaic2 


VR 


33 x 480 s 


20:39:05.8 


-18:37:45.5 


1.3 


NEP0813NW3 


Blanco 


2002-08-13 


Mosaic2 


VR 


21 x 480 s 


20:45:10 


-17:38:27.7 


0.87 




Blanco 


2002-08-15 


Mosaic2 


VR 


20 x 480 s 


20:44:58.9 


-17:39:20.7 


0.674 


NEP0815NE3 


Blanco 


2002-08-15 


Mosaic2 


VR 


21 x 480 s 


20:47:25.0 


-17:32:20.8 


0.6 



Table 4 

Observation Details. * - Field used only for confirmation purposes. A - ecliptic 
longitude (°). (5 - ecliptic latitude (°). A at - longitude with respect to Neptune (°). 



Field 


Area ( c 


'square) 


A 





Aat 


UNE 


0. 


32 


322.9 


-0.72 


15.4 


UNW 


0. 


.32 


322.9 


-0.72 


15.4 


MEGAPrimeNl 


0. 


.85 


334.9 


-0.76 


19.9 


MEGAPrimeN2 * 






334.9 


-0.76 


19.9 


N10032W3 


0. 


.38 


307.4 


0.13 


0.0 


N11033 


0. 


.38 


307.4 


0.13 


0.0 


NEP0813NW3 


0. 


.38 


309.6 


0.06 


0.0 


NEP0815NW3 * 






309.6 


0.06 


0.0 


NEP0815NE3 


0. 


.38 


309.6 


0.06 


0.0 
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Table 5 

Detection efficiency parameters using Eq. 2, and ffux measurement error parameter 
using Eq. 3. A - peak efficiency, m* - magnitude where efficiency is half the peak, g 
- width parameter. Z - field zeropoint. 7 - error parameter. 



Field 


Vr 


ncix 




g 


Z 




T 


UNE 


n 


96 


25 32 (~R) 


41 


26.21 


n 


47 


UNW 





.97 


25.44 (R) 


0.40 


26.22 


0. 


,49 


MEGAPrimeNl 





.97 


26.43 (g') 


0.41 


26.46 





.27 


N10032W3 


0. 


.91 


25.10 (VR) 


0.46 


26.0 





.77 


N11033 


0. 


.93 


25.20 (VR) 


0.5 


26.0 


0. 


.58 


NEP0813NW3 





.93 


25.13 (VR) 


0.34 


25.9 





.73 


NEP0815NE3 





.97 


25.18 (VR) 


0.27 


25.9 





.58 
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NEP0813NW3a5 *C 


23.58 ±0.06 
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43 ± 2.5 


1.4 ±0.6 


NEP0813NW3a7 *C 


24.05 ±0.10 


VR 
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VR 
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25.05 ±0.23 


VR 


42 ±3 


5±6 


NEP0815NE3a5 


25.49 ±0.36 


VR 


44 ±3 


6±6 



Table 6: Detections List. * - object used in likelihood 
fits. C - object confirmed on second night, g - object 
fell on chip gap on second night. V- object has variable 
magnitude. Number in brackets in column 2 is the range 
of magnitude measurements for the variable objects. 
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Trujillo et al. (2001) See Appendix 2 
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Trujillo et al. (2001) See Appendix 2 
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Trujillo et al. (2001) See Appendix 2 


TE5M 


g' 


7 


5 


1.44 ±0.9 
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Trujillo et al. (2001) See Appendix 2 



Table 7: Results of the maximum likelihood fits for all 
fields fit individually. Uncertainties are taken from the 
extrema of the 1-sigma likelihood contours. *-Provides a 
3-sigma lower limit only. N - Number of objects detected 
in the field. A50 - number of objects brightward of 50 % 
threshold. 
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Fig. Captions 



Fig. 1 Variation of magnitude of reference stars versus hour angle of obser- 
vations. Smooth line is the magnitude variation with airmass expected from 
nominal airmass extinction terms reported by the telescopes. The reference 
image used for planted object magnitude scaling is shown as an open triangle 
for each night, with a typical uncertainty shown for each point given. Each field 
is labeled; a: MEGAPrime Nightl, b: MEGAPrime Night2, c: N10032W3, d: 
N11033, e: UNE, f: UNW, g: NEP0813NW3 Nightl, h: NEP0815NW3 Night2, 
i: NEP0815NE3. 

Fig. 2 Grid of shift rates and angles used to search the MEGAPrime field. 
Rates are in arcseconds per hour. Angles are in degrees below the horizontal. 
The source visible in the images is a real 23.8 mag. object with a rate of motion 
close to 2 arcsec. hr _1 at 20.8°, consistent with the ecliptic. The subtraction 
well and characteristic trailing can be seen at other rates of motion. 

Fig. 3 Distribution of inclinations and distances for the artificial objects planted 
in the observations, taken from the planted objects lists of 10 chips in the 
MEGAPrime field. Solid circles mark those planted objects that were found 
during the object search. Open circles mark those that were not found. The 
same magnitude distribution of artificial objects was used at all distances. 
Thus, for a given magnitude, the detection efficiency did not depend on ob- 
ject distance. 

Fig. 4 Net efficiency for each field with 50% efficiency marked with a vertical 
line. Magnitudes are converted to R-band using average colours (see Section 6). 
Points are the binned efficiencies determined from the image search. Errorbars 
are 1-sigma Poisson confidence regions. The solid curve is the best fit efficiency 
using Eq. 2. The dotted curve is a cumulative histogram of visually rejected 
false candidates. Each field is labeled; a: MEGAPrime Nl, b: UNE, c: UNW, 
d: N10032W3, e: N11033, f: NEP0813NW3, g: NEP0815NE3. 

Fig. 5 Credible regions for the maximum likelihood fits to individual fields 
with 5 or more detections. Contours have been shifted to R-band using the 
nominal colours (see Section 6). Solid: 1-sigma, Dashed: 2-sigma, Dotted: 3- 
sigma credibility contours. Field name and number of detections used in fits 
presented. Point indicates the best fit when combining all data together. Con- 
tours indicate m is not consistent between surveys unless small shifts in m 
are allowed. 

Fig. 6 1,2, and 3-sigma confidence regions for the single power-law maximum 
likelihood fits using all data. Plotted are the likelihood contours for a and m . 

Fig. 7 Top: Histogram of combined data using 0.4 mag bin-widths. Object 
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magnitudes are shifted to R-band using typical magnitude colours. Errorbars 
are 1-sigma Poisson intervals. Solid straight line: best fit line with a = 0.65 and 
m = 23.42. Middle: Effective area (area times efficiency) of combined data. 
Dashed line: histogram of all detections using 0.4 mag bin- width. Bottom: Log 
ratio of histogram to best fit a = 0.65 and m = 23.42. Errorbars are 1-sigma 
Poisson interval. 

Fig. 8 Probability of rejection of a broken power-law from the data. Plotted 
is faint-end slope 0:2 versus break magnitude m#. Contour represents the 95% 
rejection confidence region. Contour has been smoothed to remove effects of 
course simulation sampling. 
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Fig. 2. 
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Fig. 5. 
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Fig. 7. 
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